Disentangling protostellar evolutionary stages in clustered 
environments using Spitzer-IRS spectra and comprehensive SED 

modeling 



O 
(N 

>v 



Jan Forbrich 1 , Achim Tappe 1 , Thomas Robitaille 1 , August A. Muench 1 , Paula S. Teixeira 2 , 
Elizabeth A. Lada 3 , Andrea Stolte 4 , & Charles J. Lada 1 



Oh- 
6 



> 

o 

in 

o 
o 



3 



ABSTRACT 

When studying the evolutionary stages of protostars that form in clusters, the role of any 
intracluster medium cannot be neglected. High foreground extinction can lead to situations 
where young stellar objects (YSOs) appear to be in earlier evolutionary stages than they actually 
are, particularly when using simple criteria like spectral indices. To address this issue, we have 
assembled detailed SED characterizations of a sample of 56 Spitzer-identified candidate YSOs 
in the clusters NGC 2264 and IC 348. For these, we use spectra obtained with the Infrared 
Spectrograph onboard the Spitzer Space Telescope and ancillary multi-wavelength photometry. 
The primary aim is twofold: 1) to discuss the role of spectral features, particularly those due to 
ices and silicates, in determining a YSO's evolutionary stage, and 2) to perform comprehensive 
modeling of spectral energy distributions (SEDs) enhanced by the IRS data. The SEDs consist 
of ancillary optical-to-submillimeter multi-wavelength data as well as an accurate description of 
the 9.7 /an silicate feature and of the mid-infrared continuum derived from line-free parts of 
the IRS spectra. We find that using this approach, we can distinguish genuine protostars in the 
cluster from T Tauri stars masquerading as protostars due to external foreground extinction. Our 
results underline the importance of photometric data in the far-infrared/submillimeter wavelength 
range, at sufficiently high angular resolution to more accurately classify cluster members. Such 
observations are becoming possible now with the advent of the Herschel Space Observatory. 

Subject headings: circumstellar matter — infrared: stars — open clusters and associations: individual 
(NGC 2264) — open clusters and associations: individual (IC 348) 



1. Introduction 

Finding and characterizing individual Young 
Stellar Objects (YSOs) in embedded clustered en- 
vironments is challenging because the observed 
YSO spectral energy distributions (SEDs) can eas- 
ily be influenced by dense material that is dis- 
tributed throughout the cluster, but physically un- 
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related to an individual YSO or its envelope. In 
that situation, the extra extinction can mislead 
the observer since an affected source may appear 
to be in an earlier evolutionary stage than it ac- 
tually is. For example, the SED of a T Tauri star 
behind foreground extinction can in some circum- 
stances look protostellar. In a more isolated envi- 
ronment, this situation would only occur in case 
of edge-on disks, aggravated by the presence of 
disk fl aring. To avoid confusion, iRobitaille et al 



(J2006I ) proposed that sources be assigned designa 
tions that refer to the actual evolutionary stages, 
rather than spectral classes defined by observa- 
tions, since the latter are affected by extinction. 
In most instances, the evolutionary stages (I, II, 
III) correspond to the SED classes. However, in 



heavily extincted regions the SED classes may not 
always correspond to the appropriate evolutionary 
stages. In embedded clustered environments, the 
difference between YSO classes and evolutionary 
stages becomes particularly important. 

Infrared spectroscopy allows us to go beyond 
protostellar classifications based on photome- 
try. At near-infrared wavelengths, varying de- 
grees of veiling can be used to distinguish prote- 
st ars from PMS stars dCasali fc Matthews! 1 1992 ; 
Greene fc Ladalll996l; IWhite et al.ll2007l ). At mid- 
infrared wavelengths, silicate dust and ices cause 
telltale spectral features that can be used to re- 
fine spectral classes. Many of these features oc- 
cur in the wavelength range that was accessible 
with the Infrared Spectrogr aph (IRS) onboar d 
the Spitzer Space Telescope (JHouck et al.l 120041 ). 
In the context of star formation, the power of 
this instrument has been shown most notably by 
studies of T Tauri stars and protos tars in Taurus , 
Chamaeleo n, and Ophiuchus (Fur lan et al.l 200 
20081 I2OO9I) . Most recently. iMcClure et al.l (J2oT 
study an extensive sample of YSOs in Ophiuchus, 
using IRS data to determine their evolutionary 
stages. 

To find out whether the effect of intracluster or 
external foreground extinction can be disentangled 
from the actual evolutionary stage of a YSO, we 
have obtained IRS spectra and assembled detailed 
multi-wavelength photometry for a sample of 56 
candidate YSOs in NGC 2264 and IC 348. We 
then use these data for comprehensive SED model- 
ing with the primary aim of ascertaining the evolu- 
tionary stages for the sources in our sample. These 
can then finally be compared to previous classifica- 
tion attempts that are mainly based on less broad 
mid-infrared wavelength coverage. For each of the 
56 sources, we have obtained mid-infrared spectra 
using Spitzer-IRS. These data not only allow us to 
study spectral features that trace the evolutionary 
stage of a YSO (silicates and ices), but they also 
allow us to derive photometry from line-free parts 
of the spectrum to accurately reflect the underly- 
ing continuum. This technique also allows us to 
accurately describe the main silicate spectral fea- 
ture by photometric data at selected wavelengths. 
Combined with ancillary photometric data span- 
ning optical to submillimcter wavelengths, the IRS 
data thereby allow us to construct SEDs that are 
tailored toward the requirements of detailed mod- 



eling. The accuracy of the modeling depends on 
the SEDs not containing spectral lines that are 
not part of the underlying model (e.g., PAH fea- 
tures); these features can fall into the band passes 
of broadband filters. Our IRS-dcrivcd photometry 
attempts to depict the continuum and the silicate 
features. Compared to the usual approach of us- 
ing photometry from Spitzer TRAC in SEDs for 
subsequent modeling, we discard bands 3 and 4 of 
the IRAC data that we have for every source in fa- 
vor of corresponding line-free photometry derived 
from the IRS data. 

1.1. Source selection 

We selected 56 candidate protostars in NGC 
2264 and IC348 - 40 in the former and 16 in the 
latter region - from previou s Spitzer photometric 
stu dies of these re gions by iMuench et al.l (|2007l ) 
and iTeixeiral (|2008l ). Our target sources are listed 
in Table [U 

NGC 2264 is a prominent, clustered star- 
forming region loca ted at a distance of ~913 pc 



(Bax ter et al.l 120091 ): for a recent review, see 
Dahml (|2008l ). Here, we focus on the popula- 
tion of embedded protos tars in the vicinity of 
IRS 2 (JYoung et all 120061 ) where we selected 40 
sources for Spitzer-IRS mid-infrared spectroscopic 
observations. This region is one of the two most 
prominent regions of ongoing star formation in 
NGC 2264, next to IRS 1. All selected sources 



are l ocated in the Spokes cluster (jTeixeira et al 
20Q6). Figure Q] shows the selected sources in 
a mid-infrared view of the region. IC 348 is lo- 
cated in Perseus at a distance of ^260 pc (e.g., 



Lombard i et al. | 120101); for a recen t review, see 
Herbstl (|2008l) . IMuench et al.l (|2007l) conducted a 
Spitzer census of this star-forming region. The 
IC 348 sources discussed in this paper are among 
the most embedded in the cluster and are predom- 
inantly protostars in a narrow filamentary ridge 
located about 1 pc southwest of the central clus- 
ter. In Figure [2 we mark the selected sources in 
a mid-infrared view of the region. 

A first comparison within our sample can be 
accomplished by using the IRAC spectral in- 
dex «iraCo as determined from a least-squares 
fit to the p hotom etry of all four IRAC bands 



(jLada et al.ll2006l ): see section l2~2l for a discussion 
of the photometry. The sources in our NGC 2264 
subsample show diverse IRAC spectral indices, 



ranging from a = —1.41 (source 27) to a = 4.21 
(source 45). Eighteen sources in this sub s ample 
have positive spectral indices. ISung et al.l (J2009I ) 
identify 37 out of the 40 sources as candidate 
cluster members in various stages, 20 of them are 
listed as class I sources; their classific ation is in- 
cluded in Table [TJ Note that while ISung et al. 



( 20091 ) list no membership status for the three 
remaining sources, our results suggest that they 
are cluster members, too. In the IC 348 sample, 
the IRAC spectral indices range from a = —0.93 
(source 14) to a = 1.92 (sour ce 6) . The classifi- 
cation from lLada et al.l ( 20061 ) and iMuench et al 



(|2007l ) is listed in Table [Q Most of the sources 
in this subsample are class I objects. Since we 
do not find significant differences between these 
two cluster subsamples - both (unevenly) sample 
the same YSO classes -, we will discuss them as 
one sample in the following, focusing more on a 
discussion of evolutionary stages in general. Note, 
however, that the source table allows the reader 
to identify individual YSOs in both clusters. 

In this paper, we first discuss the Spitzer- 
IRS observations and their data reduction. We 
then describe the ancillary multi-wavelength data, 
ranging from optical to submillimeter wave- 
lengths, and how we assembled composite SEDs 
from both IRS spectra and ancillary data. Sub- 
sequently, we perform comprehensive modeling of 
these SEDs before concluding with a discussion 
and summary of the results. 

2. Observations and Data Reduction 

2.1. Spectroscopic observations with Spitzer- 
IRS 

IC 348 and NGC 2264 were observed with the 
Spitzer Infrared Spectrograph (IRS) between 2006 
September 16 and 2007 October 30 (Spitzer pro- 
gram 30033). We mapped the majority of YSOs 
using the IRS short-low (SL2/SL1, 5.2-8.7/7.4- 
14.5 fim, R = 64-128), short-high (SH, 9.9- 
19.6 /mi, R = 600), and long-high (LH, 18.7- 
37.2 /mi, R = 600 ) or long-low (LL, 14.0-38.0 /mi, 
R = 64-128) modules. 

For IC 348, the sources are covered by individ- 
ual small maps of typically 2 by 3 dithered point- 
ings, i.e., three steps perpendicular and two steps 
along the slit, for both low and high resolution 
modules. We included separate high-resolution 



background observations for each source. For 
NGC 2264, which shows a much higher degree of 
source clustering, we designed extended SL/LL 
maps that cover a larger number of sources si- 
multaneously, with individual high-resolution 2 by 
3 dithered maps for the brightest sources. For 
the observations in NGC 2264, only a single high- 
resolution background observation was obtained; 
it is used for all NGC 2264 high-resolution maps. 
All maps in this program are o ptimized for 
data reduction with CUBISlVflvl.60 (jSmith et al" 
20071 ). i.e., they were obtained with half slit-width 



stepping perpendicular to the slit and additional 
dithering along the slit direction in most cases, 
resulting in quadruple data redundancy for each 
pixel in the map center. The individual steps of 
our data reduction are as follows: 

1. Acquisition of the BCDs (Basic Calibrated 
Data) as provided by the Spitzer IRS 
pipeline S17.2 to S18.2 

2. Importation of all sequential observations 
taken at the same date into CUBISM 
projects, separately for each module 

3. Selection and subtraction of an average 
background (dedicated for SH/LH; selected 
from outrigger slit positions of the same 
and/or other sources judged by variability 
of the background level in IRAC/MIPS (In- 
frared Array Camera, Multiband Imaging 
Photometer for Spitzer) maps 

4. Identification and excision of remaining 
bad pixels using CUBISM's statistical al- 
gorithms 

5. Construction of the data cube without us- 
ing CUBISM's SLCF (Slit Loss Correction 
Function), i.e., keeping the original IRS 
pipeline calibration (see fc|2.1.1j) 

6. Extraction of spectra for each source using 
a square aperture of typically 4 by 4 mod- 
ule pixels (nominal pixel sizes for the IRS 
SL/LL/SH/LH modules are 1.8, 5.1, 2.3, 
and 4.5", respectively) 

7. Conversion of spectra from the default 
MJysr -1 unit to Jy by multiplying by the 
aperture area in steradian 



1 http://ssc.spitzer.calt ech.edu/archanaly/contributed/cubism/ 



8. Impo rtation of spectra into SMART ([Higdon et al 
120041 ) for joining of the orders and scaling to 
match the MIPS 24 ^m photometry (used as 
reference, see below) 

2.1.1. CUBISM 

The Spitzer IRS pipeline is optimized for iso- 
lated stellar point sources. It is correcting for light 
losses due to the limited size of the slit and the 
extraction aperture compared to the point-spread- 
function. The flux calibration is accomplished by 
observing standard stars and matching their spec- 
tra to stellar models. Hence, the pipeline flux cal- 
ibration is invalid if: a) the observed source is not 
a stellar point source, b) the slit is not centered 
on the observed point source, or c), the default 
pipeline extraction apertures are changed. It is 
possible to correct for items b) and c) by deriving 
a new calibration using the pipeline approach. 

However, if the observed source is not a stel- 
lar or 'true' point source, e.g., a YSO with an 
accretion disk and/or an extended envelope, we 
cannot expect the flux calibration to be correct. 
To make matters worse, the effects change with 
wavelength and also among the different IRS mod- 
ules because of their different sizes and the result- 
ing uneven coverage on an asymmetric, extended 
source. Since all YSOs have an individual and 
a priori unknown extended structure and, in ad- 
dition, slit coverage will change as a function of 
distance to the source, no reliable flux calibration 
can be derived from standards. Our approach to 
this dilemma is to use CUBISM (see ij2.ip . which is 
optimized for handling Spitzer IRS mapping data 
of homogenous extended sources, and scale the re- 
sulting spectra to the measured MIPS-24 /im pho- 
tometry of each source. As a consistency check of 
the photometric calibration, we have extracteco 
color-correct magnitudes for IRAC bands 3 and 
4 from the IRS spectra in order to compare the 
results with the IRAC aperture photometry. For 
spectra that are not too noisy, the IRAC photome- 
try yields fluxes that are slightly higher than those 
from the synthetic IRS photometry, with mean 
offsets of 12% in band 3 and 18% in band 4, al- 
though offsets of up to a factor of two are occur- 
ring. While the two sets of photometry are equiv- 



2 using IDL code provided by the Spitzer Science Center at 
http://ssc.spitzer.caltech.edu/dataanalysistools/cookbook/10/ 



alent in their spectral coverage, this comparison 
does not account for slit losses. Additionally, any 
source variability would cause further differences 
since the IRAC and IRS data were not obtained 
simultaneously. Given these limitations, the two 
sets of photometry thus appear to be consistent. 

We compared the extracted spectra of CUBISM 
and SPICE for a representative, mapped YSO in 
our sample. This example, source 56, is shown in 
Figure [3J We selected this source since it is fairly 
bright and shows complex spectral structure due 
to ices, silicate grains, H2 emission lines, and a ris- 
ing continuum due to the presence of an envelope. 
SPICE is the IRS data reduction tool supported 
by the Spitzer Science Center, and it is optimized 
for single pointing (staring) observations of stellar 
point sources. The SPICE spectra were extracted 
by using only the map slits centered on the YSO, 
whereas the CUBISM spectra were extracted us- 
ing the approach described in ij2.ll Reassuringly, 
both methods yield very similar spectral features, 
continuum shapes, and signal-to noise ratios. In 
both cases additional scaling is needed to join the 
spectral orders and match the MIPS photometry. 
The main advantages of CUBISM for reducing IRS 
mapping data of a large sample of diverse YSO's 
are: 

• The data reduction is greatly facilitated due 
to integrated BCD bookkeeping, accurate vi- 
sualization and accounting of all slit point- 
ings, background subtraction, and statistical 
bad pixel handling. 

• The spatial maps cover all or a large fraction 
of an extended YSO and therefore minimize 
the effects of uneven slit coverage or varying 
slit orientation. 

• There is no restriction on the size and shape 
of the extraction aperture for the source and 
background spectra. 

• The spatial maps provide increased data re- 
dundancy compared to single slit staring ob- 
servations. 

The resulting spectra, concatenated from the 
different modules, are shown in Figures 0] to [JJ 
where they arc sorted by evolutionary stage, as 



informed by subsequent analysis! Overall, they 
represent a wide variety both in evolutionary stage 
as well an in the spectral features present. 

2.1.2. Continuum and ice features 

Many of the YSOs in our sample show promi- 
nent ice absorption features at 5-8 and 15/xm as 
well as silicate absorption at 9.7 and 18 /zm (see, 
e.g., Figure |3j). In order to quantify and com- 
pare the ice and silicate features between different 
sources, we measured the center optical depths 
t\ of the 6.0 /im and 9.7 jum H 2 ice and silicate 
absorption features, respectively, and the inte- 
grated optical depth J r dv of th e 15 ^m ice fea- 



ture which is dominated by CO2 (jGerakines et al 



19991 ) a nd has been descri b ed as a C0 2 :CO:H 2 
blend (jPontoppidan et al.l l2008h . The latter 
can be converted to a column density, A ico — 
A^l J ' t dv, where A- lce is the ice band strength of 
1.1 x 10~ 17 cm molecul e -1 assuming pure CO 2 ic e 



(JGerakines et al.l Il995c IPontoppidan et al.l |2008) 



In order to convert the spectra to an optical depth 
scale, t\ = ln(/ cont ,A^ )> we determined the con- 
tinuum by fitting a low-order polynomial on a 
log-log scale to the spectral regions below wave- 
lengths of about 5.8 /im and above 30 /J,m (see 
Figure [3} . These regions ar e least affected by ab - 
sorption features (see also iBoogert et al.l |2008) . 
In addition, we determined a local continuum for 
the 15/Ltm ice feature by fitting a low-order poly- 
nomial to the adjacent spectral continuum of that 
feature. A detailed analysis of the individual ice 
feature components is beyond the scope of this 
paper. We note that the distribution of spectral 
feature strengths is similar to previously published 
surveys of other regions, even without taking into 
account the different sets of criteria that were used 
to define the various samples. 

2.2. Ancillary optical, infrared, and sub- 
millimeter photometry 

Our IRS spectra complement a wealth of ancil- 
lary photometric data that we use for subsequent 
comprehensive S EP modeling. For IC 348, we use 



data compiled bv lMuench et al.l (|2007|) . The only 



exceptions are MIPS mid-infrared and SCUBA 



submillimeter data where we have converted sev- 
eral flux densities for individual sources into upper 
limits if there is any indication that several sources 
contribute to the flux density. 

Our main source of optical UBVRI photom- 
etry for NGC 226 4 is the anc ill ary d ata cat- 
alog compiled by JRebull et all (|2006l ). The 
photometry for the sources present ed here is 



from yarious sources c ited t h erein (|Sung et al. 
19971; iFlaccomio et al.l Il999t IPaxk et all l200q 
Rebull et al.l 120021 lLamm et al.ll2004 . To resolve 



3 Note that sources 14 and 49 are not shown because nearby 
bright sources blend with their data at longer wavelengths. 
The SEDs of these sources therefore are not modeled at all. 



inconsistencies (for a given source, the UBVRI 
photometry typically is from at least two differ- 
ent sources), the USNO-B1 catalog has been used 
additionally. 

In the ne ar-infrared, we use FLAMINGOS 
(jElstonl Il998l ) imaging data in JHK to surpass 
the 2MASS catalog in sensitivity. The observa- 
tions were taken on 01 January 2005 when the 
instrument was mounted to the 2.1m telescope 
on Kitt Peak. In this data, all but seven of our 
sources in NGC 2264 are detected at least in K 
band. For the undetected sources, the estimated 
sensitivity limits of J = 19 mag, H = 18 mag, 
and K — 17.5 mag were assumed as upper lim- 
its. Aperture photometry was performed on these 
images with an aperture of five pixels or 1.5" 
with a surrounding annulus of five pixels width 
to subtract the background. To avoid saturation 
problems, we have used 2MASS photometry for 
three bright sources (sources 17, 36, and 47). 

In the mid-infrared range covered by Spitzer, 
we have obtained aperture photometry from IRAC 
and PRF-fitting photometry from MIPS imaging 
data. The IRAC data were obtained in March 
2006 (program ID 37) when both 0.6 sec and 12 sec 
exposures were taken. Aperture photometry was 
performed with an aperture radius of five pixels 
and a background annulus with inner and outer 
radii of 5 and 10 pixels. For all but two sources, 
we used the 12 sec exposures. Sources 20 and 47 
required correcting for saturation in the first two 
IRAC channels. We used a comparison of IRAC 
and longer- wavelength data to find obvious blends 
of multiple objects and exclude the corresponding 
SEDs from modeling; this concerns sources 6+7 (a 
blend), 41, and 43. From the IRAC photometry, a 
spectral index o^rac was determined from a least- 
squares fit to all four IRAC channels in \F\. To 
ensure comparability within our full sample, the 



same method was used to also determine spectral 
indices for th e sources IC 348 , using the photomet- 



(2007). 



ric data from iMuench et al. 

In the submillimeter range, we use aper- 
ture photometry to estimate upper limits in 
SCUBA data (45 urn a nd 850 /im) obtained by 



Wolf-Chase et all (|2003l ) for NGC 2264 in 1999 
and 2000. For NGC 2264 and most of IC 348, we 
decided to only use upper limits due the crowding 
of the cluster relative to the angular resolution. 
In determining the photometry, we read out the 
pixel value at the source position (calibrated in 
Jy/beam) and subsequently added five times the 
rms of the surrounding 25 pixels. For the use in 
SED modeling, in order to use strict upper lim- 
its, we have doubled the upper limits that were 
determined in the mid-infrared (MIPS) and sub- 
millimeter bands. 

2.3. Constructing photometric SEDs 

As a preparation for the SED model fits, we 
have constructed photometric SEDs for all ob- 
served sources using the ancillary photometric 
data as well as the IRS observations. The SEDs 
cover the optical (UBVRI), the near infrared 
(JHK), the mid-infrared (3/zm up to 35 /im from 
IRAC and IRS) and the submillimeter ranges 
(450 /im and 850 /im). We have extracted pho- 
tometry from the IRS data to complement the 
broadband data from the literature, using bands 
that characterize the continuum as well as the 
9.7 /im and 18 /im silicate features - these are 
taken into account in the SED models that we 
use - while being least affected by molecular or 
atomic lines and ice features that are not taken 
into account in the model. The wavelengths for 
which photometry was derived using the IRS data 
are 5.580, 7.650, 9.95, 12.93, 17.72, 24.28, 29.95, 
and 35.06 /im, where the first two bands replace 
IRAC bands 3 and 4 since they are less affected by 
spectral features than these two IRAC bands. The 
widths of the IRS spectral windows that are con- 
verted into photometry are set to 2% of the center 
wavelengths. In this process, a few photometry 
points were flagged and not taken into account 
in the modeling if the IRS data were too noisy. 
Overall, fluxes from 21 separate wavelength bands 
were used to construct the SEDs for individual 
sources. 



3. Analysis 

3.1. IRS spectra 

3.1.1. Dust and ice features 

In the following, we will focus on the silicate and 
ice features and their correlation with the sources' 
infrared spectral indices while also noting the de- 
tection of PAH emission in only three sources in 
NGC 2264, perhaps indicating the presence of rela- 
tively massive stars. The spectral analysis for the 
sources in our sample is summarized in Table [2j 
where the sources are already sorted by evolution- 
ary stage, as informed by subsequent SED mod- 
eling. The plots of the spectra in Figures [4] to [7] 
are sorted in the same fashion. The detailed SED 
modeling that is required to sort the sources by 
evolutionary stage will be described subsequently. 

For the purpose of comparing spectral features 
to an infrared spectral index as a more easily ob- 
servable source parameter, we use the IRAC spec- 
tral index simply because it is available for more 
sources in our sample than another frequently used 
index, spanning from the near-infrared K band 
to 24 /im. Several, presumably deeply embedded 
sources in our sample remain undetected in the 
near-infrared wavelength range. We note, how- 
ever, that the trends discussed here are also ap- 
parent when using the K-24/im spectral index 
(not shown) . We remind the reader that the IRAC 
spectral index was calculated from all four IRAC 
bands; to ensure the comparability with other 
datasets, it does not use any IRS data in replace- 
ment for IRAC channels. The scatter of the four 
IRAC channels is reflected in the uncertainty of 
the spectral index, also derived from the fit. 

We now discuss the relationship of the 9.7 /im 
silicate feature and the total ice column density 
with the IRAC spectral index. We also directly 
compare silicate emission and ice absorption and 
compare this relationship to the results of the c2d 
survey. In Fig.[8j we show the corresponding plots 
in six panels that are described in more detail in 
the following. As shown in Fig.[5K/b, a source with 
a positive spectral index tends to have stronger 
silicate absorption and ice features than a source 
with a negative spectral index. Note that the sili- 
cate feature can appear in emission (indicative of 
disks) or in absorption, while the ice features al- 
ways are absorption features. Our discussion of 



silicate features focuses on the prominent 9.7 /um 
silicate feature, but we note that the optical depth 
of the 18 /xm silicate features is correlated to the 
one at 9.7 /jm, even if any more detailed compari- 
son is limited by the uncertainties. The measured 
full-width half-maximum (FWHM) line widths of 
the 9.7 /um silicate feature range from 1.8 /im to 
3.6 /im. Fig. [8h shows that the broader line widths 
correspond to negative silicate optical depths, i.e., 
silicate emission, while the line width does not 
change much for the silicate features in absorption. 
There is a tentative anti-correlation of the width 
of the silicate feature with the IRAC spectral in- 
dex; the broadest lines are found in emission and 
in the sources with the lowest spectral indices (see 
Fig.[8ji), possibly indicative of grain growth in the 



later e volu tionary stages . How ever, I Sargent et al. 



(|2009l ) and l Watson et all (|2009l ) point out that the 
effects of grain growth on the silicate feature are 
difficult to observe unambiguously even for class II 
sources where one might expect the effect to be 
more discernible than in protostars. Finally, there 
also appears to be a correlation between the op- 
tical depth of the 9.7 /Ltm silicate absorption fea- 
tures and the total ice column de nsity (shown 
i n Fig . |5b), similar to the results of 
2008, Fig. 29). As pointed out by 



Furlan et al. 



Furlan et al. 



2008), the silicate optical depth and the ice col- 



umn density are expected to be correlated since 
a higher silicate optical depth should also provide 
shielding for ice mantles to grow. While gener- 
ally, sources with larger optical depths of the sili- 
cate feature tend to have stronger ice features, the 
class I source 53 is the only source not following 
this trend; it has a much lower ice column density 
than what could be expected from this correlation, 
indicating that the central star might be able to 
evaporate the ice mantles. 

None of the sources with silicate in emission 
show ice features. This probably indicates that 
the envelopes of these sources have dissipated and 
ices are no longer observable. In Fig. |5f, we com- 
pare the silicate optical depth and the total ice 
column density from thi s study with results from 
the Spitzer c2d sur vey (jPontoppidan et al.ll2008t 
Boogert et al.l 120081 ) . Again there is a correlation 
of large silicate optical depths and high ice column 
densities even though there are two exceptions 
similar to source 53. These two sources (W33A 
and W3 IRS 5) have high silicate optical depths, 



but comparably low total ice column densities. 
W33A is an H II region and W3 IRS 5 is a massive- 
star-forming region. Possibly, both sources as well 
as our source 53 produce strong enough radiation 
to evaporate the ice mantles. 

A quantitative comparison of the spectral fea- 
tures with the literature is difficult since partic- 
ularly the methods used to determine the con- 
tinuum vary considerably. However, we can di- 
rectly compare our result s with the c2d surve y 
(jPontoppidan et al.l 120081: iBoogert et ail 120081) . 
The range and distribution of the optical depths 
of the silicate features (T6 Mm and Tg.y^m) and of 
the total ice optical depth are very similar. 

3.1.2. PAH and Gaseous Spectral Lines 

As noted above, a detailed discussion of all 
spectral features that were observed in our sam- 
ple is beyond the scope of this paper. However, we 
would like to point out the most important IRS re- 
sults beyond silicate and ice features. We detected 
strong PAH emission in sources 25, 36, and 42, all 
in NGC 2264. The emission is clearly associated 
with the YSOs, since the surrounding background 
is devoid of extended PAH emission in all three 
cases. However, the Spitzer resolution does not 
permit to resolve the origin of the PAH emission 
within the YSO envelope/disk system. We can ex- 
clude extended PAH emission beyond the Spitzer 
resolution limit for sources 25 and 42. In these 
cases, the PAH emission is hence located within a 
radius of about 1000 AU from the YSO sources. 
In the case of source 36, the high luminosity as 
well as blending of the PAH 6.2 /im feature with 
ice features make the spatial PAH analysis uncer- 
tain. These three sources might represent the up- 
per end of the mass function in the Spokes clus- 
ter. We will later see that among these sources 
only the SED of source 42 can be modeled reli- 
ably. The fit indeed suggests that source 42 is 
in the upper luminosity range among the objects 
whose SEDs we have studied in more detail. Note 
that even though NGC 2264 IRS 2 is not included 
in our spectral map, its immediate vicinity also 
shows strong PAH emission. 

Within the Spokes clu ster, lYoung et al.l (J2006I ) 
and lTeixeira et al] ([20071 ) discovered a dense clus- 
ter of protostars called the microcluster. This re- 
gion is included in our sample as two separate 
objects, 41 and 43, each being a blend of several 



sources that can be resolved in the near-infrared 
wavelength range. For these regions, we can take 
advantage of the spectral mapping capabilities of 
Spitzer to characterize the region. Even though 
the microcluster remains spatially unresolved at 
the longer wavelengths, we can still study a variety 
of spectral features in the region; this is shown in 
several panels in Fig. [9] The interpretation of the 
spatial features is complicated by the highly vari- 
able extinction across the region, most strongly af- 
fecting the near-IR wavelengths. In addition, the 
absorption maps are constrained by the presence 
of a background continuum source at the respec- 
tive wavelength, i.e., the absence of absorption 
does not necessarily indicate the absence of the 
absorbing species. 

Fig. IH1 (panel b) shows spatially extended H 2 
emission, which together with the detection of 
[Ne II] indicates outflow activity in the vicinity of 
sources 41, 43, and 44. The H2 detection is consis- 
tent with the narrowband H2 (1- 0) 2.12 /xm image 
presented by I Young et al.l ([20061 ) ■ The 30 /an con- 
tinuum emission and the 15.2 /zm ice absorption 
(Fig. El panel c) and f), respectively) are consis- 
tent with the position s of the protost a rs SM A-2, 3, 
5, and 6 reported by iTeixeira et al.l ( 2007 ). when 
taking into account the larger Spitzer beam. Note 
that SMA-5 and 6, the class protostars with the 
highest observed 230 GHz flux in the microcluster, 
have no near-IR counterpart due to the high ex- 
tinction, but are detected in our longer wavelength 
Spitzer maps. 

The aim of a quantitative analysis of spectral 
lines is usually to determine the local physical con- 
ditions, primarily density and temperature of the 
gas, and to infer chemical abundances. Infrared 
emission lines are a useful gas diagnostic for op- 
tically obscured regions such as embedded YSOs. 
However, a detailed, quantitative analysis of in- 
frared emission lines in YSOs is hampered for the 
following reasons: 

• Typically, only a limited set of infrared lines 
is observed for a given source (probing cer- 
tain ionization stages and energy levels). 

• Inadequate angular resolution results in lim- 
ited knowledge of the origin, physical geome- 
try, and location of the emissive region (out- 
flow shocks, an accretion disk, an envelope, 
etc., or a combination of those). 



• Collisionally excited lines from an inhomoge- 
neous medium will be biased toward regions 
of high temperature, where collisionally ex- 
cited emission lines will be strongest (e.g., 
shocked gas). 

• The UV radiation field (outflow shocks, disk 
and outflow cavity irradiation) and infrared 
continuum are hard to quantify, particularly 
if the location of the emitting material is un- 
known. 

• The line-of-sight extinction is difficult to 
quantify accurately. 

In light of these difficulties, a meaningful dis- 
cussion of the relative line intensities of the sources 
in our sample is beyond the scope of this paper 
and would ideally require spatially resolved ob- 
servations. However, a qualitative assessment, in 
particular of the presence or absence of infrared 
emission lines, can still give useful insights into 
the physical characteristics of a given YSO. Ta- 
bic [4] shows an overview of the infrared emission 
lines observed in our sample and lists important 
atomic and molecular line parameters. 

The last column in Table |4] gives the critical 
density, n criM = J2j<i A ij/J2j^ilij^ where i is 
the index for the initial level with energy Ei, j is 
the final level, 7 is the collisional rate coefficient 
out of level i, and A is the Einstein coefficient for 



spontaneous emission (|Osterbrocldll989. eq. 3.31). 
If the gas density is much larger than the criti- 
cal density of a given energy level i, then colli- 
sional de-excitation dominates over spontaneous 
emission, and the infrared lines originating from 
this level get collisionally suppressed. If the den- 
sity is much smaller than the critical density, the 
system exhibits a non-LTE (local thermodynamic 
equilibrium) behavior, where line emission is the 
dominant de-excitation mechanism. 

The initial energy level of a given infrared line 
is usually populated by collisions with the most 
abundant gas species (electrons, H, or H2, depend- 
ing on the conditions) and is depopulated either 
via spontaneous radiative emission giving rise to 
the observed spectral line or via collisional de- 
excitation. The general conditions necessary for 
an infrared emission line to be strong are: a) a 
large abundance of the species, b) a large popula- 
tion of the initial energy level (a kinetic gas tern- 



perature similar or larger than Ei/ks), and c) a 
large Einstein A coefficient for spontaneous emis- 
sion. Note that these are qualitative conditions, 
i.e., substantial populations and abundances can 
compensate for small Einstein A- values (e.g., H2 
pure rotational lines). 

The most common infrared lines observed in 
our sample of 56 sources are pure rotational tran- 
sitions of H2 (17 sources, of which 10 sources show 
S(4) and higher) and atomic fine-structure tran- 
sitions of Ne + (9 sources) and Fe + (7 sources). 
Among these sources, there is no clear correlation 
between the detection of any one of these features 
with evolutionary stage, partly certainly due to 
the low number of sources in each category. For 
example, out of the 17 sources with H2 transitions, 
four are in stage I, six are in stage II(ex), and three 
are in stage II (while the stages of four sources 
were not conclusively determined). 

The excited H2 lines trace the presence of 
hot molecular gas, typically from slower, non- 
dissociating outflow shocks. Good examples are 
sources 41 and 43, which are part of the micro- 
cluster, an association of young sources with clear 
signs of outflow activity (see above and Fig. [5]). 
The [Ne II] 12.8 /im line is a good example of a 
fine-structure line with a relatively large A-value, 
a low- lying energy level, and a large critical den- 
sity. In addition, neon is not depleted in the gas 
phase since it is not incorporated into dust grains. 
Therefore, Ne + is abundant if sufficiently ener- 
getic UV photons are present to ionize neutral 
neon (A < 57 nm). This may be the case in the ir- 
radiated surface layers of accretion disks or in fast 
outflow shocks that are strong enough to ionize 
atomic hydrogen. 

Iron is strongly depleted in the gas phase and 
therefore it is typically observed in fast, dissocia- 
tive shocks, which destroy dust grains and liber- 
ate the iron into the gas phase (e.g., SNR shocks 
and fast protostellar outflows). Fe + has a com- 
plex energy level structure and features both low 
and high-lying energy levels, which makes it a very 
useful density and temperature diagnostic. 

Combining all the information above, we at- 
tempt a simplified interpretation scheme: 

H2 0—0 S(4) and higher traces a hot molecular 
environment (elevated gas temperature of a 
few 1000 K or larger, low UV radiation, non- 



dissociating shock). 

[Ne II] 12.8 /.tm traces a high UV radiation en- 
vironment leading to high ionization and H 2 
photodissociation. 

[Fe II] lines trace faster, dissociating shocks that 
partially destroy dust grains and produce 
UV radiation/ionization. 

Of the eleven sources in our sample that show 
[Ne II] and/or [Fe II] lines, five show lines of both 
species, four show only [Ne II] and two show only 
[Fell]. The [Nell] 12.8 /xm emission in sources 
with no observed [Fe II] lines is unlikely to be due 
to strong outflow shocks. These are probably cases 
of irradiated disks with little or no outflow shock 
activity, consistent with the stage II/II(ex) clas- 
sification of three out of these sources (while we 
did not determine a stage for the fourth source). 
The two sources with [Fe II] lines but no [Ne II] 
12.8 /im both seem to be young stage I sources, 
one of them being part of the microcluster, with 
strong outflow shock activity. 

Sources that show both [Ne II] and [Fe II] lines 
are probably due to irradiated disks and/or strong 
outflow shocks. For example, source 41, which is 
part of the microcluster and thus has no assigned 
stage due to source crowding, shows strong H2, 
[Fe II], and [Ne II]. It is unlikely that H 2 and Ne+ 
coexist, and indeed our spectral mapping indicates 
a spatial offset between the [Nell] 12.8 /im and 
the extended H2 emission, consistent with an ir- 
radiated disk and an outflow (see Fig.[SJ panel b). 
Note as a caveat that the near-infrared imaging 
in the figure shows the complex, blended emission 
in the microcluster and the clear limitations of the 
spatial resolution of our mid-infrared observations. 
An advantage is, however, that the longer wave- 
length observations penetrate the highly obscured 
region and show the deeply embedded protostellar 



3.2. SED Modeling 

3.2.1. Procedure 



The SEDs of all the sources were fit 



the la rge set of model SEDs from iRobitaille et al. 



using 



(120061 hereafter R06) and the SED fitting tool 



from IRobitaille et all ([2007, hereafter R07). These 



models were computed using the radiation trans- 



fer code from IWhitnev et al.l (|2003bl lah. The ba- 
sic model consists of a pre-main-sequence star sur- 
rounded by a flared accretion disk and a rotation- 
ally flattened envelope with cavities carved out by 
a bipolar outflow. Both scattering and reprocess- 
ing of the stellar radiation by dust, as well as the 
disk accretion luminosity, are included. The set 
of models covers a large range of stellar masses 
(from 0.1 to 50 M ) and evolutionary stages (from 
deeply embedded protostars to 'anemic' and 'tran- 
sitional' disks). In total, 14 parameters describe 
the stellar properties and the circumstellar dust 
geometry, such as the stellar luminosity, temper- 
ature, disk mass, radius, inner radius, envelope 
infall rate, or bipolar cavity opening angle. For 
each set of parameters, SEDs were computed at 
10 different viewing angles. In total, 20,000 sets 
of parameters, and therefore 200,000 model SEDs 
were computed. 

The R07 SED fitting tool fits each of the mod- 
els in the R06 model set to the data, allowing both 
the distance and external foreground Ay to be free 
parameters (within specified ranges) and charac- 
terizes each fit by a x 2 value. By choosing a x 2 
threshold for 'good' fits, one can then look at how 
well constrained different parameters are. 

For the fitting, we ignored models with view- 
ing angles i > 80° (i.e., edge-on). Model fits 
with edge-on models tend to require unrealistically 
high stellar luminosities and seriously affect the re- 
sulting parameter ranges for many of the sources. 
However, if the sources in our sample are randomly 
oriented, we would expect 56/9 ~ 6 sources to 
have i > 80°. However, this number is an up- 
per limit because edge-on sources are fainter, and 
therefore a flux-limited sample such as ours will 
be biased towards non-edge-on sources. There- 
fore, while we may be incorrectly fitting a couple 
of edge-on sources with non-edge-on models, we 
prevent the majority of the sources from being fit 
by edge-on models, which seriously skew the re- 
sulting parameter ranges. 

In this paper, we define the goodness-of-fit cri- 
terion to be x 2 ~ Xbost < 3 x "data- This criterion 
is identical to that chosen in R07, and is arbitrary, 
but models that fit better than this criterion gen- 
erally fit the data well. While this criterion is fairly 
'loose' in statistical terms, as described in R07, 
the coarse coverage of 14-dimensional parameter 
space, intrinsic uncertainties in the models, and 



variability all mean that adopting a stricter cri- 
terion would lead to over-interpretation. For the 
purposes of the modeling, we allow the external 
Ay to vary between and 40 mag and the distance 
to vary within 5% ranges of 913 pc for NGC2264 
and 260 pc for IC 348. 

The interstellar extinction from the sun to the 
outer 'edge' of the star forming regions is not likely 
to exceed a magnitude of visual extinction. As will 
be shown further on in this paper, the external 
foreground (intracluster) extinction inside both 
star forming regions is what dominates the overall 
extinction toward most of the sources. Therefore, 
rather than assuming the same dust model as in 
R07, we choose to use a dust model with larger 
grains (M. Wolff, private communication, 2009). 
The dust model contains a mixture of astronomi- 
cal silicates and graphite i n solar abundance, us - 
ing the optical constants of lLaor fc Draind (|1993l ) , 
and a size distrib ution derived from a maximum 
entro py method (]Kim et al.l Il994t IClavton et al.l 



2003) to reproduce an Ry value of 5.5. Such high 



Rv values are more suited to dense star form- 
ing regions tha n the typical diffuse ISM value of 
R v =3-3.5 (e.g.. ICardelli et al.|[l989h . In order to 
assess the influence of the extinction law that we 
assume for the external foreground extinction on 
the source classification, we have also used the two 
observational mid-infrared extinction laws derived 
bv lMcClurel ([2009). As noted in Table El we find 
that only six to seven out of 44 sources would be 
classified differently. There seems to be no system- 
atic effect on the number counts of the three differ- 
ent evolutionary stages since the changes largely 
cancel out. Conclusions concerning source pop- 
ulations are thus not significantly affected. The 
choice of extinction law does, however, alter the 
parameters derived for individual YSOs from the 
SED fits, but these are not considered in detail in 
this paper. 

Figures ITU1 to IT21 show the SEDs of most of the 
sources along with the model fits. These sources 
are very well fit by the models, and the derived pa- 
rameters are reliable (within the framework of the 
assumed model geometry). On the other hand, for 
the sources shown in Figure IT3l the derived param- 
eters cannot be trusted. In the case of sources 8 
and 32, strong ice features in the 10 to 20 /im range 
and the lack of longer wavelength emission means 
that the model fits cannot be trusted. Sources 
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12, 16, 25, and 36 are not well fit by the models. 
In particular, source 25 is likely to be an edge-on 
source. This source, which has strong PAH fea- 
tures, has a very similar overall SED and mid- 
infrared spectrum as Gomez' Hamburger which 
is a k nown edge- on dis k (e.g., iBuiarrabal et al. 



20081 : IWoodet alJl2008ft . Finally, in the case of 
sources 6+7 (a blend), 41, and 43, the broadband 
fluxes at most wavelengths contain contributions 
from multiple sources. Therefore, while the SEDs 
appear well fit, the resulting parameters are not 
meaningful. 

Table [3] lists the lower and upper values of each 
parameter of the good fits for each source. The 
best-fit value is not meaningful, and is therefore 
not shown. Only sources with reliable parameters 
(as discussed above) are listed. The modeling re- 
sults are discussed in the next section. We note 
that the fit results are generally compatible with 
additional constraints from the literature (spectral 
types, X-ray-derived column densities); deviations 
are not significant. 

3.2.2. Results 

One particularity about the R06 set of mod- 
els is that parameter space was not sampled in a 
uniform way due to computational limitations. In- 
stead, the parameter space coverage was designed 
to include reasonable regions of parameter space 
suggested by observations and theory. Therefore, 
there exist intrinsic correlations between the pa- 
rameters in the set of models that need to be 
taken into account when usi ng the models to look 
for trends (|Robitailld 120081 ). For example, if in 
this sample we were to find an inverse correlation 
between stellar temperature and envelope mass - 
which would make sense as one expects the stel- 
lar temperature to rise and the envelope mass to 
decrease as the source evolves - such a trend was 
assumed in the R06 set of models and therefore 
the result may just reflect this assumption rather 
than a real trend. We find no correlation between 
model parameters that is not likely to be due to 
the assumptions built into the R06 models and 
therefore in the remainder of this section we con- 
centrate on the determination of the evolutionary 
stage of the sources. 

The sources in the present sample were orig- 
inally targeted as preferentially being candidate 
embedded protostars, based in part on the spectral 



indices derived from Spitzer IRAC and MIPS pho- 
tometry. As shown in R06, while such a spectral 
index provides a way of statistically estimating the 
evolutionary stage of sources in a sample, it is not 
necessarily reliable on a source-by-source basis es- 
pecially in regions that can be expected to contain 
large amounts of intersource extinction such as 
massive, cluster-forming cloud cores. Therefore, 
some of the sources in this embedded cluster sam- 
ple may not actually be embedded protostars. In 
this section, we use the results of the SED model- 
ing along with information extracted from the IRS 
spectra and the spatial distribution of the sources 
to examine in detail the physical nature of each of 
the sources. 

Based on the results of the broadband SED 
modeling, we separate the sources into sources 
that require significant envelopes, and are there- 
fore are likely to be genuine protostars, and ones 
that can be fit by disk-only models, and therefore 
may be more evolved. We refer to the former as 
stage I sources and the latter as stage II sources 
as in R06, but define the separation between the 
two groups of sources differently to R06. We de- 
fine the genuine embedded sources (stage I) as 
those which cannot be fit by models which have 
envelope masses less than 10 _3 M Q , and the po- 
tentially more evolved sources (stage II) as those 
that can be fit with models that have envelope 
masses below 10 _3 M Q . Stage I sources are ex- 
pected to display class or class I SEDs. Typ- 
ically, we would expect stage II sources to dis- 
play class II SEDs, but in the presence of signif- 
icant (unrelated) foreground extinction they can 
display class I-like SEDs. The same effect can be 
caused by flared disks that are seen edge-on, but 
this should only affect a minor fraction of a given 
sample. 

As can be seen in Table [3j some sources require 
very little external foreground extinction while 
others require at least 10 or 20 magnitudes of 
such extinction in addition to the self-extinction 
by the circumstellar dust. We therefore separate 
the stage II sources into ones that require an ex- 
ternal extinction of at least 7 mag (stage II(ex)), 
and sources that do not require an external Ay of 
more than 7 mag (stage II). Stage II(ex) sources 
typically display flat-spectrum/class I-like SEDs 
while the stage II sources typically display more 
traditional class II SEDs. The exact values to sep- 
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arate the sources into the three categories sources 
is fairly arbitrary, but there are only a few bor- 
derline cases. Most sources clearly fall into one of 
these categories. We note that some stage II(ex) 
sources have a minimum interstellar extinction of 
less than 7 mag: these are cases where only the 
stage I models require external extinctions below 
7 mag, and all stage II models require higher ex- 
ternal extinctions. 

Figure Q3] shows the spatial distribution of the 
sources overlaid on a SCUBA 850 /im image for 
each cluster. One striking result in NGC 2264 
is that the stage II sources clearly do not follow 
the same spatial distribution as the other sources, 
and are preferentially located to the West. The 
stage I sources are preferentially associated with 
brighter SCUBA 850 /im emission, while the stage 
II sources are preferentially associated with unde- 
tected or faint SCUBA 850 /im emission, with the 
stage II(ex) sources associated with intermediate 
brightness sub-mm emission. The mean 850 /im 
fluxes associated with the sources and the cor- 
responding standard errors of these means are 
0.69±0.10, 0.34±0.06, and 0.14±0.06 Jy/beam for 
stages I, II(ex), and II respectively. 

Taken together, these results suggest that the 
stage II sources are a population of T-Tauri-likc 
stars that are in their protoplanetary disk phase, 
and furthermore are not deeply embedded in the 
molecular cloud, as they show fairly low extinc- 
tions. On the other hand, the sources in stage 
II (ex) require at least 7 mag of external foreground 
extinction and are spatially coincident with the 
massive dense SCUBA core in the East. This indi- 
cates that the sources in the East are more deeply 
embedded in the cloud than those in the West. 
Both the fact that the stage II (ex) sources can be 
fit by disk-dominated models and their lower as- 
sociated SCUBA 850 /im fluxes suggest that these 
are in fact a similar population to stage II, but 
with larger amounts of external foreground extinc- 
tion. Finally, the stage I sources are genuine em- 
bedded protostars. 

4. Discussion 

After the analysis of the IRS spectra and sub- 
sequent detailed SED modeling, we can character- 
ize the YSOs in our sample more accurately than 
what was previously possible from YSO spectral 



classes alone by assigning evolutionary stages next 
to their SED classes. As we have seen, the compre- 
hensive SED modeling suggests that the sources 
fall into three classes: 1) genuine protostars (stage 
I), 2) disk sources behind high foreground extinc- 
tion (stage II(ex)), and 3) disk sources with less 
external foreground extinction (stage II). Of the 
44 sources for which we determine evolutionary 
stages, 17 are stage I protostars. Most of these 
were already classified as class I in the literature 
(one, source 4, is a class source, but a stage ear- 
lier than stage I has not been defined), with the 
exceptions of sources 28, 45, and 51 which have not 
been classified, or have be e n clas sified differently 
(source 35) in ISung et al.l (J2009I ). On the other 
hand, as summarized in Table [I] the sample con- 
tains many previous class I sources that instead 
appear to be in stage II(ex). Thus, a number of 
sources appear to be in later evolutionary stages 
than what was previously suspected. Most previ- 
ously determined class II sources turn out to be in 
stage II, although a few are in stage II(ex). 

With regard to the use of spectral indices as an 
indicator of YSO class, our sample contains a few 
objects that can serve as instructive examples of 
the limits of using a single spectral index for the 
classification of individual sources. Sources 3 and 
56 both have negative IRAC spectral indices, but 
it turns out that they are stage I protostars and 
were classified as class I protostars in the litera- 
ture. 

The IRS data are a significant help to constrain 
the SEDs of these sources both by allowing us to 
make use of spectral properties in YSO charac- 
terization, and by extending the SED wavelength 
coverage. The latter is further improved since the 
IRS data allow us to avoid spectral ranges contam- 
inated by unwanted line features for photometry 
and subsequent SED modeling. In particular, the 
youngest sources show strong silicate and ice fea- 
tures in their mid-infrared spectra; both features 
are correlated with the sources' spectral indices, 
both IRAC and K-24/im. A correlation of the 
silicate optical depth at 9.7 /im and the total ice 
column density is present, but has more scatter 
than the correlations of the two features with the 
IRAC spectral index as a measure of evolutionary 
stage. The main exception appear to be objects 
with high silicate optical depths but low total ice 
column densities. These might be more massive 
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stars with radiation fields able to evaporate the ice 
mantles. A similar correlation i s found in the c2d 
data and bv lFurlan et al.l (|2008l ). We detect three 
sources with PAH emission (all in NGC 2264), a 
small subset of our sample. Presumably, these are 
more massive YSOs. For the crowded region of 
the microcluster of NGC 2264, we can make use 
of IRS spectral maps, for example to identify re- 
gions with jet activity. 

Based on the IRS and modeling results, we can 
check whether the three evolutionary stages that 
were identified by SED modeling in turn have dis- 
tinctive features in their IRS spectra. We have 
therefore marked the sources in Fig. |8] according 
to their evolutionary stage. Indeed, the stage I 
sources tend to have large optical depths of the sil- 
icate feature and strong ice features. Stage II(ex) 
sources have less strong silicate and ice features 
and stage II sources tend to show the silicate fea- 
ture in emission while not showing any ice fea- 
tures. It is interesting to note that the IRS spec- 
tral features also appear to be distinctly differ- 
ent when comparing stage II (ex) and stage II 
sources. Fig.[8]shows that while virtually all stage 
II sources have the 9.7 jam in emission, indica- 
tive of disks, that is not the case for the stage 
II(ex) sources, presumably due to the additional 
extinction which may 'cancel out' the emission. 
Also, none of the stage II sources show ice fea- 
tures, while stage II(ex) interestingly do show ice 
absorption, although at smaller total column den- 
sities than the stage I sources. In our interpreta- 
tion, the stage II(ex) ice features would be related 
to the additional extinction by intracluster fore- 
ground material and not the sources themselves. 
In summary, while IRS data are of crucial impor- 
tance, they alone are not sufficient to distinguish 
the evolutionary stages. For this purpose, detailed 
SED modeling is required. 

We point out that in the context of our compre- 
hensive SED modeling, the archival submillime- 
ter data were very important for the separation of 
the sample into the different evolutionary stages. 
Dense dust and gas in an envelope manifest them- 
selves in submillimeter maps as concentrated emis- 
sion in contrast to more large-scale foreground ex- 
tinction. Future high-angular resolution observa- 
tions in the far-infrared/submillimeter wavelength 
range, as forthcoming from the Herschel Space 
Observatory, will be even better suited to help de- 



termining the evolutionary stages in YSO samples. 
Finally, we discuss the observable quantities 
that correlate best with the different evolutionary 
stages. We have seen that this is certainly not a 
spectral index since the correlation of the evolu- 
tionary stages with spectral indices is not unam- 
biguous. The genuine protostars (stage I) all have 
submillimeter detections and the highest silicate 
optical depths as well as the highest ice column 
densities. The stage II sources are broadly char- 
acterized by a silicate feature in emission and no 
ice features. The stage II(ex) sources lie between 
these two extremes, but their SEDs can be eas- 
ily mistaken for class I sources (in fact, some were 
previously determined to be class I sources, see 
above); in particular, they do not show the sili- 
cate feature in emission. Once an individual SED 
shape is well characterized, far-infrared and sub- 
millimeter continuum data are the most important 
photometric data to distinguish the different evo- 
lutionary stages. 

5. Summary and Lessons Learned 

We have combined Spitzer IRS spectra, SCUBA 
submillimeter, Spitzer mid-infrared, FLAMIN- 
GOS near-infrared and published UBVRI pho- 
tometry with detailed radiative transfer models 
to investigate the nature of the most embedded 
members in the NGC 2264 and IC348 clusters. 
We summarize our results as follows: 

• Empirical SED classes and mid-infrared 
spectral features do not provide unam- 
biguous identifications of YSO evolutionary 
stages in clustered environments. 

• In particular, significant amounts of fore- 
ground extinction can transform a class II 
SED into a flat/class I SED and suppress 
silicate emission characteristic of disks. 

• SED modeling incorporating the silicate fea- 
tures and mm/submm photometry provides 
a more reliable measure of the evolutionary 
stage. 

• In our sample of IC 348 sources, we identify 
three stage I, six stage II (ex), and one stage 
II source. In our NGC 2264 sample, covering 
the Spokes cluster, we identify 14 stage I, 9 
stage II(ex), and 11 stage II sources. Among 
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the previously known 20 class I sources in 
our NGC 2264 sample, we find 16 stage I 
sources whereas the reminder are classified 
as stage II(cx). 

• Only in NGC 2264, we identify three sources 
with PAH emission, potentially tracing the 
upper end of the mass function in this clus- 
ter. 

• We find various YSOs with rotational H2 
transitions as well as sources with neon and 
iron fine structure transitions, but for these 
features, we do not find a clear correlation 
with evolutionary stage in our sample. 
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Table 


1 














Source 


LIST 










ID 


RA 


Dec 


oirac 


C«K-24 


class a 
(lit.) 


Stage 
(this paper) 


1 


LRL 00245 


03:43:45.17 


+32:03:58.7 


0.44 


0.42 


I 


II(ex) 


2 


LRL 54362 


03:43:50.97 


+32:03:24.0 


0.95 


0.04 


I 




3 


LRL 54361 


03:43:51.03 


+32:03:07.7 


-0.83 


n.a. 


I 




4 


LRL 57025 


03:43:56.89 


+32:03:03.4 


n.a. 


n.a. 







5 


LRL 00013 


03:43:59.64 


+32:01:54.2 


-0.31 


-0.27 


thick 


II(ex) 


6 


LRL 54459 


03:44:02.41 


+32:02:04.5 


1.92 


n.a. 


I 


b 


7 


LRL 54460 


03:44:02.62 


+32:01:59.6 


0.97 


n.a. 


I 


b 


8 


LRL 01916 


03:44:05.78 


+32:00:28.5 


-0.29 


0.00 


I 




9 


LRL 04011 


03:44:06.91 


+32:01:55.4 


-0.15 


0.20 


I 


II(ex) 


10 


LRL 00051 


03:44:12.97 


+32:01:35.4 


0.37 


0.42 


thick 


II(cx) 


11 


LRL 01889 


03:44:21.35 


+31:59:32.7 


-0.01 


0.58 


I 


II(ex) 


12 


LRL 00435 


03:44:30.28 


+32:11:35.3 


0.47 


0.51 


Hal 




13 


LRL 01872 


03:44:43.31 


+32:01:31.6 


0.40 


0.76 


I 


II(ex) 


14 


LRL 01898 


03:44:43.89 


+32:01:37.4 


-0.93 


0.65 


0/1 


C 


15 


LRL 00234 


03:44:45.22 


+32:01:20.0 


-0.42 


-0.32 


I 


II 


16 


LRL 00904 


03:45:13.81 


+32:12:10.1 


1.20 


0.53 


I 




17 


SSB 08807 


06:40:39.34 


+09:34:45.6 


0.03 


0.01 


II.+ 


II 


18 


SSB 08910 


06:40:40.06 


+09:35:02.9 


-0.66 


-0.58 


H,+ 


II 


19 


SSB 09077 


06:40:41.15 


+09:33:57.9 


-1.02 


-0.63 


II.+ 


II 


20 


SSB 10207 


06:40:48.62 


+09:35:57.7 


-0.04 


0.13 


I.+ 


II(ex) 


21 


SSB 10281 


06:40:49.15 


+09:37:36.4 


-0.57 


-0.24 


pre-TD 


II 


22 


SSB 10387 


06:40:49.89 


+09:36:49.5 


-1.31 


-0.89 


H,+ 


II 


23 


SSB 10455 


06:40:50.29 


+09:34:53.7 


2.63 


2.45 






24 


SSB 10587 


06:40:51.06 


+09:35:22.2 


1.73 


1.30 






25 


SSB 10710 


06:40:51.91 


+09:37:55.9 


1.89 


0.92 






26 


SSB 10769 


06:40:52.34 


+09:34:57.7 


1.75 


1.60 






27 


SSB 10775 


06:40:52.38 


+09:34:31.4 


-1.41 


0.41 




II(ex) 


28 


SSB 11017 


06:40:54.19 


+09:36:46.2 


1.45 


1.44 


member? 




29 


SSB 11314 


06:40:56.17 


+09:36:31.0 


-0.48 


-0.42 


n,+ 


II 


30 


SSB 11347 


06:40:56.40 


+09:35:53.3 


-1.16 


-0.42 


ii.+ 


II 


31 


SSB 11445 


06:40:57.00 


+09:33:01.4 


0.52 


-0.04 


ii, p 


II 


32 


SSB 11592 


06:40:57.98 


+09:36:39.5 


1.20 


1.62 


I,X 




33 


SSB 11598 


06:40:58.01 


+09:36:14.6 


-0.45 


0.44 


l,x 


II(ex) 


34 


SSB 11802 


06:40:59.14 


+09:36:14.9 


0.83 


1.27 


I 


II(cx) 


35 


SSB 11825 


06:40:59.27 


+09:33:25.1 


-0.48 


0.07 


pre-TD 




36 


SSB 11829 


06:40:59.30 


+09:35:52.4 


-0.94 


-0.61 


II, X 




37 


SSB 11847 


06:40:59.37 


+09:33:33.5 


-1.22 


-0.02 


IL + 


II 


38 


SSB 12070 


06:41:00.98 


+09:32:44.5 


-1.21 


-0.89 


IL + 


II 


39 


SSB 12208 


06:41:01.83 


+09:34:34.3 


-0.35 


0.53 


I,x 


II(ex) 


40 


SSB 12583 


06:41:04.25 


+09:34:59.7 


-0.06 


0.77 


I 




41 


SSB 12780 


06:41:05.56 


+09:34:08.0 


0.21 


1.64 


I.X 




42 


SSB 12820 


06:41:05.77 


+09:35:29.6 


0.62 


0.36 


I,X 




43 


SSB 12875 


06:41:06.19 


+09:34:08.8 


0.33 


1.53 


I,x 




44 


SSB 12893 


06:41:06.29 


+09:33:50.0 


1.09 


1.38 


I 




45 


SSB 12913 


06:41:06.42 


+09:35:54.5 


4.21 


2.58 


member? 




46 


SSB 12951 


06:41:06.66 


+09:33:57.8 


-0.08 


0.42 


I,x 


II(ex) 


47 


SSB 12966 


06:41:06.74 


+09:34:45.9 


-0.38 


0.05 


IL + 


II 


48 


SSB 12968 


06:41:06.77 


+09:33:34.9 


2.29 


2.56 


I 




49 


SSB 13069 


06:41:07.40 


+09:34:54.9 


-0.98 


0.20 


II, X 


C 


50 


SSB 13111 


06:41:07.67 


+09:34:19.1 


1.63 


2.24 


I,X 




51 


SSB 13238 


06:41:08.61 


+09:35:42.7 


1.78 


2.41 


member? 




52 


SSB 13244 


06:41:08.63 


+09:36:03.4 


-0.10 


0.35 


II 


II(ex) 


53 


SSB 13427 


06:41:09.90 


+09:35:40.8 


2.08 


2.76 


I 




54 


SSB 13499 


06:41:10.43 


+09:34:18.6 


-0.35 


0.06 


II, X 


II(ex) 


55 


SSB 13567 


06:41:10.93 


+09:34:08.3 


-0.58 


0.15 


II,x 


II(ex) 


56 


SSB 13696 


06:41:11.84 


+09:35:31.5 


-0.47 


0.43 


I,X 






NGC 2264 fr 




1 120091) with the followin 


g symbols: 


+ — Halphe 




a l'o 


omlSung et al 


emission star 



1G 



with X-ray emission; X — X-ray emission star; x — X-ray emission candidate; P — X-ray emission 
star w ith strong Hct emission; for IC 348 from |Lada et a i. (2006, disk types) and iMuench et al.| 
I J2007I) 

Sources 6 and 7 form a single blended source at most wavelengths, and arc thus discussed as 
such. 

c Due to blending with nearby sources, insufficient photometry at longer wavelengths to warrant 
modeling and further discussion. 
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Table 2 
Spectral features and fit results 



Source Stage tqo 



FWHM 9.7 D 

[cm ] 



JVtot(ice) b 
[xl0 18 cm~ ; 



Comments 



2 




1.01* 


2.79* 




0.65 








3 




1.21* 


3.00* 




2.08 


H 2 S(l)-S(7), [Nell] 12.8, 


[Fe II] 17.9, 26.0 




4 










<3.64 


H 2 S(l)-S(2) , [Fe II] 17.9 


24.5, 26.0, 35.3, [S I] 25.2, [Si II] 34.8; 


no SL 


23 




0.49 


1.85 


2.2 


2.02 








24 




1.61* 


1.94* 






no LL/LH 






26 




0.93 


2.86 


2.1 


3.25 








28 




<0.75 


1.89* 


2.1 




no LL/LH 






35 




<0.14 


<0.05 




<0.15 








40 




0.31 


1.36 


2.0 


0.85 








42 










0.78 


PAH 






44 




0.52 


1.60 


2.1 


1.41 








45 




1.41 


5.93 


2.0 


5.27 








48 




1.34 


3.08 


2.1 


2.65 








50 




0.92 


3.87 


2.1 


4.46 








51 




<1.26 


<1.58 




<0.42 


H 2 S(l)-S(2) 






53 




1.94 


6.27 


1.8 


1.72 








56 




0.68* 


1.28* 


2.1 


1.82 


H 2 S(2), S(5), [Nell] 12.8 


[Fe II] 17.9, 26.0, [Si II] 34.8 




1 


II(ex 


<0.07 


-0.35 


3.0 


<0.13 








5 


II(ex 


0.16 


0.43 


2.6 


0.23 








9 


II(ex 


<0.88 


<0.50 




<2.66 


H 2 S(l)-S(5) 






10 


II(ex 


0.13 


1.48 


2.3 


0.30 








11 


II(ex 


0.22* 


0.83 


2.8 


0.75 








13 


II(ex 


0.20 


0.45 


2.1 


0.81 








20 


II(ex 


<0.05 


<0.06 




<0.11 


H 2 S(l)-S(2), [Nell] 12.8 






27 


II(ex 


<1.08 


<0.32 






no LL/LH 






33 


II(ex 


<0.42 


<0.54 






H 2 S(2)-S(7); no LL/LH 






34 


II(ex 


0.35 


1.27 


2.3 


0.56 


H 2 S(l)-S(3), [Nell] 12.8, 


[Fe II] 17.9, 26.0 




39 


II(ex 


0.19* 


0.38* 


2.3 


0.27* 








46 


II(ex 


0.22* 


0.68* 


1.9 










52 


II(ex 


0.51* 


1.13* 


2.2 




no LL/LH 






54 


II(ex 


0.23* 


0.51* 


2.5 


0.88* 


H 2 S(l)-S(7) 






55 


II(cx 


<0.07 


<0.50 




<0.33 


H 2 S(0)-S(5), [Sill] 34.8 






15 




<0.20 


0.54 




<0.54 








17 




<0.07 


-0.25 


3.6 


<0.17 


H 2 S(0)-S(3), [Nell] 12.8 






18 




<0.08 


-1.09 


2.8 




no LL/LH 






19 




<0.21 


<0.38 






no LL/LH 






21 




<0.11 


-0.54 


2.8 


<0.63 


H 2 S(0)-S(2), [Nell] 12.8 






22 




<0.06 


-0.85 


2.7 




H 2 S(2)-S(3); no LL/LH 






29 




<0.08 


-0.69 


3.0 


<0.18 








30 




<0.18 


<0.14 




<0.54 








31 




<0.02 


-0.53 


3.0 


<0.25 








37 




<0.32 


-1.03 


2.4 




no LL/LH 






38 




<0.08 


-0.60 


3.4 


<0.26 








47 




<0.08 


-0.34 


2.8 


<0.09 









H 2 S(l)-S(6), [Nell] 12.8 

H 2 S(l)-S(7), [Nell] 12.8, [Ne III] 15.5, [Fell] 17.9; no LL/LH 
PAH 

PAH 

H 2 S(l)-S(7), [Nell] 12.8, [Fell] 5.3, 17.9, 26.0; part of microcluster 
H 2 S(3)-S(7), [Fell] 5.3, 17.9; part of microcluster 
no LL/LH 



6+7 

8 

12 

14 

16 

25 

32 

36 

41 

43 

49 



0.71 


1.96 


2.4 


1.59 


0.23 


1.00 


2.6 


1.36 


<0.78 


<0.23 




<0.42 


1.11* 


2.25* 




<1.00 


<0.58 


0.82 




<0.59 


<0.68 






<0.19 


0.47 


1.18* 


2.1 


1.75 


<0.12 


<0.26 




0.29 


0.60 


1.60 


1.9 


1.72 


0.97 


2.66 


2.4 


2.19 


0.19 


0.77 


2.1 
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Note. — Silicate optical depths and ice column densities are generally estimated to be accurate to ~ 25% except for some cases, marked by 
asterisks, where the error is estimated to be ~ 50%. 

a (. . . ) — no silicate band or no data for FWHM fit; for ices: (...) — no spectral data (e.g., SL or SH/LH only) 

h (cotwo : CO : H 2 0) 

c fi2 transitions are H2(0— 0) transitions. 
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Table 3 
Ranges of parameters from the modeling with the R07 SED fitting tool. 



X-ray A,, 
(mag) 



< L ©) 



(K) 



Sp. Type 



(Mq) 



Jw infall 
(M /yr) 



<M ) 



Stage a 



to 

o 



18 
19 



36.8 
20.8 

0.0 
37.0 

7.3 

0.2 
34.2 

0.9 
23.7 
10.6 
21.8 
20.0 
28.5 
17.0 

0.0 
27.0 
10.8 



40.0 
39.7 
40.0 
40.0 
40.0 

0.2 
40.0 

3.3 
36.9 
23.0 
38.3 
40.0 
40.0 
40.0 
40.0 
40.0 
27.5 



52.9 
50.2 

481.9 
55.9 
43.7 
77.4 
50.0 
4.5 
30.8 
19.0 
28.6 

124.0 
77.3 
68.8 

476.3 

153.4 
27.2 



78.5 
52.0 
700.2 
65.4 
462.6 
77.4 
70.8 
15.2 
37.1 
24.5 
13.1 
192.3 
104.4 
620.6 
2383.3 
164.2 
37.8 



16.85± 1.45 



50.50± 5.75 



0.1 

0.2 
1.2 
1.6 
0. 1 
1.3 
0.4 
5.9 
20.6 
21.2 
2.9 
3.1 
4.9 
3.7 
8.9 
18.7 
16.6 



0.5 
0.2 
2.5 
2.3 
3.8 
1.3 
0.9 
16.8 
117.5 
114.1 
19.1 
6.9 
15.8 
15.5 
24.1 
34.6 
22.3 



2807 
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The fit results listed here, including the stages listed here were determined using an cxth 
the following differences: For the low-A K law, sources 23, 31, 40, and 44 become stage II(ex 



ccs should be treated with 1 



as the silicate 



law by M. Wolff, as d 
1 sources 10 and 27 be 



feature is not fit by the : 



ibed in the text. Using the two extinction laws derived by[McClurc | 2009:) leads to 
e stage I; when using the high-Aj^- law, almost the same changes occur, but source 



NOTE. — For each parameter, the minimum and mai 

were determined from spectral fits to NGC 2264 so 
types are from SIMBAD. 



um value are those of the models that provide a good fit (x — xt 
es by the ANCHORS project (cxc. harvard. edu/ANCHORS) : 
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Table 4 
Data for gaseous infrared transitions observed in this work 



Line ID 



Avac 

[ Mm 



IE" 

[eV] 



EJki 
PCI 



An 

Is" 1 ] 



[cm 



H 2 0-0 S(0) 
H 2 0-0 S(l) 
H 2 0-0 S(2) 
H 2 0-0 S(3) 
H 2 0-0 S(4) 
H 2 0-0 S(5) 
H 2 0-0 S(6) 
H 2 0-0 S(7) 

[Fe II] 



[No II] 
[No III] 

[Si] 

[Si II] 



28.22 
17.03 
12.28 
9.66 
8.03 
6.91 
6.11 
5.51 

25.99 
35.35 
5.34 
17.94 
24.52 

12.81 

15.56 

25.25 
34.82 



7.9 
7.9 
7.9 
7.9 
7.9 

21.6 
41.0 



510 
1015 
1682 
2504 
3474 
4586 
5830 
7196 

554 
961 
2694 
3496 
4083 

1123 
925 

570 

413 



2.9E-11 
4.8E-10 
2.8E-09 
9.8E-09 
2.6E-08 
5.9E-08 
1.1E-07 
2.0E-07 

2.1E-03 
1.6E-03 
6.4E-05 
5.9E-03 
3.9E-03 

8.6E-03 
6.0E-03 

1.4E-03 

2.1E-04 



1 

13 
77 
279 
766 
1781 
3614 
6612 

7665 
4274 
297 
16185 
10133 

311319 
113670 

93131 

784 



"ionization energy needed to reach the respective ionization 
stage 

critical d ensity of initial energy level i: using H— H 2 collision 
rates for H 2 (Wrathmall ct al. 2003), electron collision rates for 
atom ic lines (Fe : Ramsbottorn etal.| [20071 . Nell: I Griffin et a! 
20011. Nelll: Bell 2001, SI: Tayal 20Ql Sill: IDufton fc Kingston 
1994J) , and a temperature of 2000 K 
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Fig. 1. — Objects selected in NGC 2264 for IRS survey, overlaid on a MIPS 24 fim image. The boxes mark 
the location of spectral maps in the indicated modules. The large spectral map covers the Spokes cluster. 
The brightest source is NGC 2264 IRS 2, and it is not included in our sample. 
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Fig. 2. — Objects selected in IC348 for IRS survey, overlaid on a MIPS 24 /im image. 
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1000 500 

Wavenumber (cm) 



Fig. 3. — Top: Spitzer IRS low and high resolution 
spectrum of source 56: CUBISM (color) compared 
to SPICE (grey/black) and IRAC/MIPS photom- 
etry (colored circles). No scalings have been ap- 
plied. Center: Continuum fit for Spitzer spectrum 
of source 26. The IRS low- and high-resolution 
modules have been scaled and merged. Vertical 
dashed lines mark the regions used for the con- 
tinuum fit. Bottom: Optical depth spectrum for 
source 26. 
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Fig. 4. — IRS spectra of the 17 identified genuine protostars (stage I sources), in \F\, with orders of 
magnitude indicated (in ergs^ 1 cm -2 ). 
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Fig. 4. — continued 
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Fig. 5. — IRS spectra of the 15 identified disk sources with high external foreground extinction (stage II(ex) 
sources) . 
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Fig. 6. — IRS spectra of the 12 identified disk sources with low external foreground extinction (stage II 
sources) . 
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Fig. 7. — IRS spectra of sources where the modeling was found to be unreliable, for example due to blending 
of several sources (see text). 
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Fig. 8. — a) Correlation of the IRAC spectral index and the optical depth of the 9.7 /^m silicate feature. 
Triangles indicate upper limits, b) Correlation of the IRAC spectral index and the total ice column density. 
Triangles indicate upper limits, c) Correlation of the optical depth and the FWHM lincwidth of the 9.7 /im 
silicate feature, d) Correlation of the IRAC spectral index and the FWHM linewidth of the 9.7 /xm silicate 
feature, e) Correlation of the optical depth of the 9.7 /an silicate feature and the total ice column density, 
without upper limits, f) Correlation of the optical depth of the 9.7 ^m silicate feature a nd the total ice 
column density in the Spitzer c2d survey (jPontoppidan et all l2008t iBoogert et alj I2008I) , without upper 
limits. In all color plots, red denotes stage I sources, green denotes stage II(ex) sources, blue denotes stage 
II, and black denotes sources where no evolutionary stage was determined. 
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Fig. 9. — Spectral maps for the microcluster region. Wavelength regions: a) Ks 2-2.3 /um (includes H2 (1-0) 
S(0)-S(2) and H I Br 7), b) H 2 (0-0) S(7)+S(5)+S(3)+S(2), with [Ne II] 12.8 ^m contour overlaid (5a), c) 
30 minus 22 fim continuum emission, d) 6.0+6.7 /jm ice absorption (H2O and NH4"), e) silicate 9.7 /xm minus 
adjacent continuum, f) 15.2 /im ice absorption (C02:CO:H20). Boxes denote the original pixel size, i.e., 
about the actual angular resolution, fo r the Spitzer data. In panel c, the submillimeter continuum is marked 
(5a contour from iTeixeira et al.ll2006t note that only sources 41 and 43 were covered in that observation). 
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Fig. 10. — Observed SEDs of the 17 genuine protostars (stage I sources), with model SEDs providing a good 
fit overplotted. The black circles show the observed fluxes, the triangles show upper limits, the black line 
is the best fit model SED from R06, and the grey lines are other good model SED fits from R06 (where a 
good fit is defined as one with \ 2 — Xbest < 3 x ridataj- The dashed line shows the contribution of the stellar 
photosphere, including the effect of foreground extinction. Note that the data points shown reflect the input 
for the fitting; upper limits therefore are not necessarily due to non-detections (see text). 
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Fig. 10. — continued 
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Fig. 11. — Observed SEDs of the 15 disk sources with high external foreground extinction (stage II(cx) 
sources) , with model SEDs providing a good fit overplotted. The symbols and lines are as in Figure [lOj 
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Fig. 12. — Observed SEDs of the 12 disk sources with low external foreground extinction (stage II sources), 
with model SEDs providing a good fit overplotted. The symbols and lines are as in Figure 1101 
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Fig. 13. — Observed SEDs of sources where the modeling results cannot be trusted, as described in the 
text. Best fit results are only shown for guidance; note, for example, that good fits were obtained for the 
unresolved microcluster sources 41 and 43. The symbols and lines are as in Figure [TU] 
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Fig. 14. — Position of the sources overplotted on 850 /im SCUBA maps of NGC 2264 (upper panel) and 
IC 348 (lower panel). The symbols denote stage I (circles), stage II (ex) (squares), and stage II (triangles) 
sources respectively. Grayscale steps are in steps of 0.2 Jybeam -1 up to 2 Jybeam -1 . 
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